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SUMMARY 

This  report  describes  an  application  of  Clustering  and  Differential  Geometry  to  the 
characterization  of  digital  terrain  maps  by  different  surface  types  and  critical  regions.  Of 
particular  interest  is  the  determination  of  specific  terrain  types,  where  they  occur,  and 
describing  them  in  ways  which  are  invariant  to  viewing  position.  The  development  of  robust 
and  efficient  segmentation  and  associated  region  characterization  procedures  is  reported. 
Good  results  have  been  obtained  for  characterizing  digital  terrain  maps.  This  work  supports 
the  development  of  techniques  for  generating  optimal  flight  paths. 
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INTRODUCTION 


Digital  terrain  elevation  map  data  are  comprised  of  elevation  (range)  data  for  a  rectilinear 
grid.  Each  grid  element  comprises  one  arc  second  of  longitude  and  latitude  on  the  Earth's 
surface.  Methods  for  mission  planning  using  these  data  are  computationally  intensive.  If 
terrain  maps  can  be  categorised  into  higher-level  features,  such  as  hills,  valleys  and 
planar  regions,  then  the  search  for  an  optimal  trajectory  through  the  terrain  would  be 
greatly  reduced.  It  was  felt  that  the  interpretive  techniques  of  machine  vision  should  be 
applicable  to  digital  terrain  maps  as  they  are  essentially  pixel  image  data,  where  elevation 
is  equivalent  to  a  grey  scale  or  range  information. 

This  document  describes  the  results  of  an  investigation  of  the  utility  of  Computational 
Vision  techniques  for  the  characterization  of  digital  terrain  information,  the  automatic 
reduction  of  terrain  maps  to  (multi-scaled)  region  types  and  relational,  or  topographical, 
graphs.  The  specific  objective  was  the  automatic  description  of  regions  in  terms  of  what 
and  where  they  are  in  terms  which  are  invariant  to  viewing  position.  The  correct 
labelling  of  these  regions  is  critical  for  optimal  flight  planning  for  military'  operations. 

The  w'ork  described  here  falls  into  two  areas; 

•  the  investigation  of  robust  and  efficient  segmentation  and  region  labelling 
procedures,  and 

•  region  feature  extraction  and  generation  of  symbolic  descriptions  in  terms  of 
relational  graphs. 

Over  the  past  decade  several  techniques  have  been  developed  for  the  segmentation  and 
labelling  of  surfaces  for  3D  object  recognition.  These  techniques  are  based  upon 
Differential  Geometry  [1,2],  Clustering  [2]  and  the  integration  and  extension  of  current 
Evidence-Based  approaches  to  surface  and  object  recognition  [3,4].  These  techniques 
include  a  detailed  set  of  surface  feature  extractors,  range  segmenters,  rule  generation 
procedures  and  matching  algorithms. 

Many  of  these  techniques  have  been  implemented  as  part  of  the  Image  and  Pattern 
Recognition  Scheme  (IPRS).  IPRS  is  a  library  of  signal  processing,  image  processing, 
pattern  analysis,  pattern  recognition  and  object  recognition  functions.  IPRS  can  be 
broken  down  into  functional  areas  that  overlap  to  some  degree.  The  system  is 
documented  in  [5].  It  currently  comprises: 

•  FM  -  File  and  memory  routines:  saving  and  loading  images,  allocating  and  freeing 
memory, 

•  IP  -  Image  processing:  including  format  conversions,  filtering,  feature  extraction, 
compression  and  other  early  vision  processes, 

•  CL  -  Data  clustering  techniques;  standard  and  new  clustering  procedures  relevant 
to  pattern  recognition  and  rule  generation,  and 


-1- 


•  OPR  -  Object  and  Pattern  Recognition  procedures;  techniques  developed  for  the 
recognition  of  structures  invariant  to  their  rigid  motions. 

In  the  interpretation  of  a  map  by  a  mission  planner,  significant  features  of  navigational 
importance  are  identified  before  the  process  of  logistic  constraint-based  planning.  The 
identification  of  features  in  the  map  corresponds  to  the  labelling  problem  in  machine 
vision.  This  report  describes  the  adaptation  of  various  computer  vision  techniques  to  the 
specific  problem  of  labelling  regions  of  digital  terrain  range  data  which  correspond  to 
imponant  topographical  features  for  trajectory  generation. 


2  REGION  LABELLING 

Terrain  data  give  the  elevation  of  a  point  on  a  grid,  with  respect  to  mean  sea  level, 
usually  in  a  longitudinal  and  latitudinal  format.  Thus,  elevation  data  are  essentially  the 
complement  of  range  data  in  the  terminology  of  machine  vision  and  object  recognition, 
for  navigation  purposes,  the  data  are  view-dependent  insofar  as  the  depth  information  is 
determined  as  a  function  of  the  sensor  position.  Altitude  above  the  terrain  is  a  function  of 
where  you  are  in  the  terrain  grid  and  how  far  above  it  you  are  flying. 

This  type  of  surface  representation  is  called  a  Ivlonge  patch  where  depth(z)  is  a  function 
of  the  projection  plane  coordinates  (x,y) : 

z=  f(x,y). 


For  flight  path  planning,  it  is  most  important  to  have  a  representation  for  the  shape  of 
region  parts  which  is  invariant  to  view  direction  -  given  that  the  region  is  visible. 
Fortunately,  the  surface  Mean(H)  aind  Gaussian(K)  curvatures  have  such  characteristics, 
as  has  been  well  known  for  over  a  century.  Not  only  are  they  invariant  to 
parameterizations  of  a  surface  (in  this  case,  viewer  position),  but  they  are  also  invariant  to 
rigid  motions  of  the  surface  itself  [6].  However,  there  are  several  methods  of  calculating 
these  curvatures  and  some  are  more  reliable  and  representative  than  others  [Ij. 

The  general  strategy  used  to  estimate  curvatures  involves  fitting  a  surface  to  a 
neighbourhood  or  w'indow  centred  about  each  pixel,  determining  the  first  and  second 
order  partial  derivatures  of  the  surface  and  then  using  these  derivatives  in  equations  for 
the  Mean  (H)  and  Gaussian  (K)  curvatures  [7].  We  have  used  the  quadratic  surface 
fitting  procedure,  which  is  the  most  common  one  used  in  the  literature  and  which  has  the 
advantage  of  ease  of  computation  of  H  and  K  for  every  NxN  window  [1].  This  procedure 
is  described  in  Appendix  A. 

Table  1  shows  the  eight  different  region  types  determined  as  a  function  of  the  signs  of  the 
H  and  K  shape  descriptors.  The  region  type  label  is  assigned  to  the  pixel  in  the  middle  of 
the  window. 
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K>0 

K=0 

K<0 

H<0 

peak 

ridge 

saddle 

H=0 

(none) 

flat 

minimal  surface 

H>0 

pit 

valley 

saddle  valley 

Table  1:  FipHt  Fundamental  Surface  Types  from  Curvature  Signs 


3  SEGMENTATION 

Segmentation  is  a  form  of  clustering  where  pixels  are  grouped  together  in  ways  which 
reflect  the  inherent  structures  of  the  data.  The  determination  of  such  structures  can  often 
depend  upon  the  application.  In  map  interpretation  this  involves  determining  whether  a 
pixel  is  part  of  a  feature  of  interest  such  as  a  valley  or  hill.  For  instance,  consider  Figure 
1  which  shows  part  of  a  terrain  map.  It  is  fairly  obvious  to  a  human  observer  that  the 
middle  section  is  a  valley,  where  the  darker  pixels  represent  lower  elevation.  In  this  case, 
the  aim  of  segmentation  is  to  assign  each  pixel  within  a  group  to  a  label  representing  a 
valley. 

The  issue  of  segmentation  for  range  data  interpretation  has  received  a  good  deal  of 
attention  in  recent  years  [1,2].  Common  to  most  approaches  is  the  development  of 
surface  part  clustering  in  terms  of  similarities  in  surface  point  position  normals  (ie  (x,y,z) 
coordinates),  curvature  information  or  surface  curve  fitting  parameters.  Segmentation, 
in  these  low-level  terms,  does  not  guarantee  the  derivation  of  pans  which  are  consistent 
with,  for  example,  topographical  regions  defined  by  other  processes.  There  have  been 
some  attempts  to  split  and  merge  such  initially  segmented  regions  to  be  consistent  with 
known  patch  boundary  feature  bounds;  for  example,  object  parts  in  the  database  [4].  The 
nature  of  our  problem  led  us  to  explore  different  segmentation  techniques  which  vary  as  a 
function  of  the  degree  to  which  positionfwAcre)  and/or  shape]  w/wr)  information  is  used. 

For  navigational  purposes,  the  eight  features  extracted  from  the  H  and  Ks  are  too  fine. 
All  the  pilot  requires  are  the  following  three  general  features: 


•  planar, 

•  valleys  and 

•  ridges. 

The  eight  fundamental  features  found  in  Table  1  were  grouped  to  reflect  these  broad 
features;  seeTable  2. 
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K>0 

K=0 

K<0 

H<0 

ridge 

ndge 

planar 

H=0 

(none) 

planar 

planar 

H>0 

valley 

valley 

valley 

Table  2:  Three  Broad  Features  revelant  to  Navigational  Purposes  determined  from 

the  Curvature  Signs 


3.1  Segmentation  by  Shape 

Segmentation  on  shape-alone  (what),  ie  calculating  the  H  and  K  values  for  every 
window,  does  not  produce  realistic  or  useful  parts  for  this  type  of  problem.  This  is  panly 
because  shape  does  not  consider  distance  information,  per  se.  Also,  the  data,  being  so 
variable,  require  significant  smoothing  before  computations  of  curvatures  -  which,  in 
turn,  transform  the  original  data  into  a  terrain  envelope.  Figure  1  shows  the  unprocessed 
terrain  map  which  is  used  for  testing.  Figure  2  shows  the  results  of  labelling  the  regions 
according  to  the  mean  and  Gaussian  curvatures  listed  in  Table  1  for  window'  sizes  3.  5 
and  13.  Figure  3,  on  the  other  hand,  present  the  results  when  labelling  the  image 
according  to  Table  2,  w'hich  contain  more  useful  regions  for  the  task  of  labelling  features 
in  digital  terrain  maps.  However,  neither  result  is  satisfactory'.  Segmentation  on  shape- 
alone  is  not  sufficient  for  this  problem.  The  similarity  calculated  simply  in  terms  of 
curvature  values  does  not  guarantee  the  necessary  contiguity  in  positional  information  to 
obtain  reliable  labelling  of  features. 

3.2  Segmentation  by  Position 

For  segmentation  based  on  position-alone,  a  segmentation  procedure  based  on  clustering 
was  explored  (see  acknowledgements).  In  particular,  agglomeration  was  used  where  each 
surface  point  is  defined  by  (x,y,az)  coordinates,  where  a  corresponds  to  a  range 
sensitivity  factor.  Agglomeration  groups  sample  points  according  their  distances  in 
feature  space,  in  this  case,  their  scaled  coordinates  [8.9).  Again,  the  icciiaiquc  is 
recursive  in  so  far  as  points  are  joined  in  terms  of  their  relative  distances;  the  data  set 
being  ultimately  reduced  to  one  point  in  the  following  way.  First,  the  nearest  two  points 
are  fused  into  one  centroid,  from  this  new  reduced  data  set,  the  process  is  repeated  -  until 
there  is  exactly  one  point  remaining.  The  associated  tree  defines  clusters  at  different 
levels.  From  a  large  number  of  runs  it  was  empirically  determined  that  the  most  useful 
value  for  a  was  about  20.  However,  a  major  concern  with  agglomeration  is  that  it 
computationally  expensive  and  therefore  is  not  a  viable  solution  for  large  terrain  maps. 
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3.3  Segmentation  bv  Position  Using  Slicing 


In  order  to  improve  the  labelling  of  images  and  to  place  a  greater  emphasis  on  height,  the 
terrain  map  can  be  sliced  according  to  elevation.  This  makes  sense  according  to  the 
physical  interpretation  of  the  terrain  in  the  context  of  aircraft  navigation.  Here,  the  map 
is  sliced  according  to  height  by  determining,  for  example,  equally  sized  intervals 
spanning  the  whole  range.  Each  slice  or  interval  reflects  the  amount  of  terrain  within  a 
particular  interval  of  height.  The  slicing  intervals  can  be  adapted  to  fit  the  range 
resolution  required.  Each  contiguous  area  within  each  slice  is  treated  as  a  separate  entity. 
Subsequently,  when  these  slices  are  re  combined  into  a  single  map.  discrete  contiguous 
areas  remain  separate. 

The  steps  used  for  combining  the  calculation  of  H  and  K  and  slicing  are: 

•  calculate  the  H  and  K  values  for  each  N*N  window  in  the  image  and  assign 
them  to  the  mid-point  of  the  window, 

•  slice  the  terrain  map  into  X  slices, 

•  uniquely  label  the  contiguous  regions  within  each  slice. 

•  recombine  slices  into  one  map  retaining  unique  region  labels. 

•  calculate  the  average  H  and  K  values  for  each  uniquely  labelled  contiguous  area, 
and 

•  assign  the  map  features  to  unique  areas. 


Though  the  solution  presented  here  is  not  optimized,  we  have  found  that  equal  slices  of 
height  produce  results  which  are  consistent  with  the  regions  perceived  to  be  present  in  the 
data.  This  method  has  also  the  advantage  ot  having  the  ability  to  control  the  number  of 
regions  via  the  slice  resolution.  However,  instead  of  equal  slices  of  depth,  slices  could  be 
determined  using  some  heuristic  based  on  the  physical  interpretation  of  the  data. 

Figure  4  shows  results  of  slicing  the  map  into  3  slices.  Figure  5  shows  the  resulting 
recombined  map  where  each  uniquely  labelled  contiguous  area  is  coloured  differently. 
Finally,  Figure  6  shows  the  map  labelled  according  to  the  labels  found  in  Table  2  and 
determined  from  the  average  curvature  values  within  each  region  for  a  window  size  of  3. 
5  and  13,  respectively. 

The  slicing  procedure  gave  similiar  results  to  the  use  of  (x.y,20z)  in  the  clustering 
procedure  as  described  in  Section  3.1,  where  a  scaling  factor  was  necessary  to  emphasise 
height.  The  magnification  of  z  demonstrates  the  importance  of  the  actual  range  value  in 
the  formation  of  region.  However,  the  computation  requirements  for  determining  the 
slices,  labeling  each  contiguous  region  within  a  slice  uniquely  and  then  recombining  the 
slices  into  one  map,  was  found  to  be  minimal  compared  to  agglomeration. 

Figure  6  bears  the  closest  resemblance  to  the  original  terrain  map.  With  this  result,  many 
co-located  regions  (as  shown  in  Figure  5)  are  of  the  same  type  and  after  labelling,  the 
unique  areas  appear  to  be  in  the  same  region  (as  seen  in  Figure  6).  Overall  structures  and 
regions  can  be  clearly  identified  and  located  using  this  combined  slicing  and  feature 
extraction  method. 
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DISCUSSION 


At  this  stage  we  have  shown  how  one  elementary  descnption  in  terms  ot  connected 
regions  can  be  uiained  from  the  data.  However,  funher  work  needs  to  be  done  to  bring 
these  procc.  .^es  to  completion.  In  panicular.  we  need  to; 

•  improve  the  position-only  clustering  procedures. 

•  develop  range-merging  procedures  for  merging  different  regions  at  different 
heights, 

•  improve  the  region  labels, 

•  check  the  ability  of  our  system  to  compare  labels  across  view  s.  and 

•  investigate  the  tv-pes  of  pan  features  necessary  to  define  regions  and  their 
relations. 

The  next  imponant  step  to  be  taken  with  this  work  is  to  develop  more  advanced  svmbolic 
labelling  procedures  to  fit  in  with  the  path  planning  problem.  W'e  would  like  to  produce 
skeleton  paths  through  the  various  regions.  These  paths,  representing  the  various  regions, 
could  then  be  used  to  determine  the  optimal  path  through  a  terrain  map.  By  using  these 
paths,  not  every  pixel  within  an  image  needs  to  be  searched  or  examined,  and  thus, 
computational  costs  will  be  greatly  reduced. 

The  resulting  skeletons  representing  the  various  regions  could  then  be  used  in 
determining  the  optimal  flight  path  across  the  terrain. 

One  solution  to  generating  these  paths  involves  formulating  the  problem  in  state-space 
[10].  Each  state  is  described  by  the  aircraffs  position,  orientation  and  speed,  while  the 
cost  associated  with  the  state  is  a  function  of  threat,  terrain  and  the  aircraft's  manoeuvre 
required  to  change  the  state.  The  resulting  graph  can  then  be  searched  using  A*  search. 
Thus,  if  the  terrain  can  be  categonsed  in^o  higher-level  features,  such  as  hills,  valleys  and 
planar  regions,  then  the  search  space  for  an  optimal  trajectorv-  through  the  terrain  would 
be  greatly  reduced. 

Funher  work  in  sketelonisation  is  required  to  exploit  the  techniques  in  the  automatic 
interpretation  of  maps  in  mis.sion  planning.  That  is.  to  exploit  graph  search,  labelled 
regions  need  to  be  reduced  to  labelled  lines  and  venices.  This  work  has  commenced. 


5  CONCLUSIONS 

In  this  study  we  have  determined  the  feasibility  of  producing,  from  input  terrain  maps,  a 
symbolic  topographical  description  of  the  data.  A  digital  terrain  map  can  be  regarded  as 
an  image  made  up  of  pixels.  Thus,  a  number  of  computer  vision  techniques  can  be  used 
to  identify  terrain  features  for  navigational  purposes. 

We  have  found  that  the  most  efficient  and  reliable  method  involves  combining  the 
calculation  of  the  mean  (H)  and  Gaussian  (K)  curvatures  and  slicing.  Segmenting  the 
map  with  respect  to  positional  information  and  then  labelling  each  contiguous  region  for 
surface  types  results  in  the  what  and  where  topographical  map.  We  have  shown  that 
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machine  vision  techniques  have  application  on  data  types  outside  their  onginal 
conception  in  object  recognition. 
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Figure  !,  Original  terrain  map  028*128  pixels) 


Figure  2.  Map  labelled  according  to  H  and  K  listed  in  Table  1,  for  windows  3^  &  13,  respectively 


KPirnsr. 


Figure  3.  Map  labelled  according  to  H  and  K  listed  in  Table  2,  for  windows  3^  &  13,  respectively 


Figure  4.  Map  sliced  into  3  slices 


Figure  5.  Slices  recombined  into  one  map,  with  contiguous  regions  labelled  uniquely 


Figure  6.  Labelled  recombined  map  according  to  average  H  and  K,  according  to  Table  2 


APPENDIX  A 


Calculating  the  Mean  (H)  and  Gaussian  (K)  Curvatures 

The  approach  used  in  this  quadratic  fitting  procedure  is  the  following: 

•  given  discrete  sample  data,  a  continuous  differentiable  function  that  'best'  fits  the 
data  is  to  be  determined,  and 

•  compute  the  derivatives  of  the  continuous  function  analytically  and  then  evaluate 
them  at  the  corresponding  discrete  points. 

Ideally,  we  would  like  to  fit  one  smooth  surface  to  all  the  data.  This  is  not 
computationally  feasible.  Instead,  a  local  surface  fit  is  determined  for  each  NxN  window 
in  the  map.  A  local  quadratic  surface  model  is  used  for  this  purpose.  Experimental 
results  have  shown  that  this  approach  is  adequate'. 

For  a  given  window  size  (scale)  the  (least  squares)  best  fitting  quadratic  surface  range 
value  is  computed,  over  an  NxN  window  at  position  (x,y),  as: 


f'(x,y)  =  Ijj  ajj  0i(x)  ({)j(y) 


where  ajj  are  selected  to  minimize  the  error  term: 


e  =  I^xy  (  ^W)  -  f(x,y)  )^- 


The  following  discrete  onhogonal  polynomials,  (Jis,  provide  the  quadratic  surface  fit: 


0(u)  =  1 


4>i(u)  =  u 


<t)2(u)=  {  u2-(M(M+1))/3  ) 

for  M  =  (N-l)/2,  ie  half  the  window  size.  The  ajj  terms  are  determined  by: 

ay  =  Zxyf(x,y)bi(x)bj(y) 


'  Bcsl,  P.  &  Jain,  A.,  1986,  Invariant  surface  characteristics  of  3D  object  recognition  in 
range  images.  Computer  Vision,  Graphics  and  Image  Processing,  33,  33-80 
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where  the  b  terms  are  normalized  versions  of  the  onhogonal  polynomials  (1 ): 


bo(u)  =  1  /  N 

bi(u)  =  3/  (M(M+1)(2M+1)U  ) 
b2(u)  =  1/  P(M)  X  {  u2-  (M(M+1)  /  3)  } 


where 


P(M)  =  8/45  m5  +  4/9  +  2/9  -  1/9  -  1/15  M  . 


Thus,  computing  the  a  terms  using  odd  sized  windows  is  simple,  because  the  bs  are 
precomputed  for  any  window  size.  From  the  a  terms,  the  first  and  second  partial 
derivative  estimates  are  then  given  by: 


fx  =  HO, 
fy  =  aoi 


fxy  =  an. 


fxx  =  2  a20,  and 

fyy  =  2  ao2. 


The  residual  squared-error  image  is  determined  by: 


e(x,y)  =  (  f'(x,y)  -  f(x,y)  )2. 


It  is  important  to  note  that  such  smoothing  results  in  errors  at  discontinuities  or  region 
boundaries.  However,  as  we  are  concerned  with  broad  region  categories  and  their 
location,  this  has  no  significant  effect  in  the  computations  of  surface  shape. 
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The  H  and  K  values  can  then  be  determined  using  the  following  equations; 


(fxx  fyy  ^xx  fy“  "‘‘^yy^x^  '  2  fx  f y  ^xy) 

_ and 

(l+fx2  +fy2)3/2 

.  -f  2 

XX  ‘yy  ^xy 


(1  +  fx2  +  fy2)2 


These  equations  apply  to  a  Monge  patch  (view-dependent  depth  map)  case  where 
refers  to  panial  differentiation  of  f  with  respect  to  u  and  v  and  f(x,y)  to  the  view- 
dependent  range  image.  Such  computations  are  enacted  after  the  initial  quadratic  surface 
fitting  procedure. 
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